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DEAR EDITOR, 


The composition and diversity of the human vaginal microbial 
community have been investigated intensively due to the 
diversity-stability relationship (DSR) -based hypothesis for 
bacterial vaginosis (BV) etiology, which was first proposed in 
the 1990s and has received renewed interest in recent years. 
Nevertheless, diversity changes (scaling) across individuals in 
a cohort or population have not yet been addressed, which is 
significant both theoretically and practically. Theoretically, 
biodiversity scaling is the core of biogeography, and 
practically, inter-subject heterogeneity is critical for 
understanding the etiology and epidemiology of human 
microbiome-associated diseases such as BV. Here we applied 
the diversity-area relationship (DAR), a recent extension to 
the classic species-area relationship (SAR), to study diversity 
scaling of the vaginal microbiome by reanalyzing reported 
data collected from 1107 postpartum women. The model 
used here characterized the power-law (or its extension) 
relationships between accrued diversity and areas (numbers 
of individuals), upon which four biogeographic profiles were 
thus defined. Specifically, we established the DAR profile 
(relationship between diversity scaling parameter and so- 
termed diversity order (q)), similarly pair-wise diversity overlap 
(PDO) profile, maximal accrual diversity (MAD) profile, and 
ratio of individual-level to population-level diversity (RIP) 
profile. These four profiles offer valuable tools to assess and 
predict diversity scaling (changes) in the human vaginal 
microbiome across individuals, as well as to understand the 
dynamics of vaginal microbiomes in healthy women. 

The human vaginal microbiome is a complex ecosystem 
that plays critical roles in maintaining host health. As the first 
defense of the reproductive tract, the vaginal microbiome is 
critical for the prevention of opportunistic pathogen 
colonization and viral infection. For example, endogenous, 
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healthy vaginal microbiota can help protect against HIV 
infection by activating local and systemic inflammation; 
however, microbiota associated with BV can also increase 
susceptibility to HIV infection (Buvé et al., 2014; Petrova et al., 
2013). For pregnant women, the vaginal microbiota is not only 
associated with maternal health but also that of neonates, with 
the composition of the newly colonized microbiome playing a 
key role in newborn immunity and metabolic development 
(Cox et al., 2014; Dominguez-Bello et al., 2010; Olszak et al., 
2012; Rutayisire et al., 2016). Furthermore, babies delivered 
by cesarean section can have a higher risk of metabolic and 
immune diseases than those delivered vaginally (Dominguez- 
Bello et al., 2010; Sevelsted et al., 2015; Younes et al., 2018), 
although Chu et al. (2017) noted that delivery mode does not 
influence microbiome composition in newborns. Moreover, in 
pregnancy, vaginal dysbiosis is hypothesized to be a 
contributor to spontaneous preterm birth (Freitas et al., 2018; 
Romero et al., 2014a; Stout et al., 2017) and miscarriage 
(Ralph et al., 1999). 

In many healthy women, the vaginal microbiota is 
dominated by Lactobacillus spp. (Macklaim et al., 2013; Ravel 
et al., 2011). Several studies (Brotman et al., 2014; Gajer et 
al., 2012; Ma & Li, 2017; Ravel et al., 2011) have confirmed 
the five major community state types of the vaginal 
microbiome in adult women, as first identified by Ravel et al. 
(2011). Four types are dominated by Lactobacillus spp., 
including L. iners, L. crispatus, L. gasseri, and L. jensenii. 
However, 20%-30% of asymptomatic, otherwise healthy 
individuals lack lactic acid bacteria in their vaginal 
microbiome, which instead consists of obligate anaerobic 
bacteria (Ravel et al., 2011, 2013). In addition, the frequency 
of microbiome type varies in different ethnic groups, with 
those microbiome not dominated by Lactobacillus spp. more 
commonly found in healthy Hispanic and black women than in 
Asian or white women (Ma et al., 2012; Ravel et al., 2011). 
Furthermore, the composition of the vaginal microbiome is 
dynamic during life and associated with menopause stage 
(Muhleisen & Herbst-Karlovetz, 2016). Recent research 
demonstrated the vaginal microbiome of perimenarcheal 
adolescents to be dominated by Lactobacillus spp., including 
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L. crispatus, L. iners, L. gasseri, and L. jensenii, similar to that 
found in reproductive-age women (Hickey et al., 2015). In 
premenopausal women, the vaginal microbiota is still 
dominated by L. crispatus and L. iners, but Lactobacillus spp. 
are often replaced by Streptococcus and Prevotella in the 
perimenopausal and postmenopausal stages (Brotman et al., 
2014). Shifts in vaginal microbiome have also been observed 
during and after pregnancy. For example, diversity and 
richness of the vaginal microbiome is lower in pregnant 
women than in non-pregnant women (Freitas et al., 2017). 
Furthermore, Romero et al. (2014b) showed that the vaginal 
microbiome of pregnant women contains a higher abundance 
of L. vaginalis, L. crispatus, L. gasseri, and L. jensenii, and a 
lower probability of switching to a Lactobacillus-deficient 
community. In addition, radical changes in Lactobacillus-poor 
vaginal communities have been found at delivery, which can 
persist for up to a year (DiGuilio et al., 2015). 

Despite extensive studies on the human vaginal 
microbiome, what constitutes normal or healthy vaginal 
microbiota remains unresolved. For example, Doyle et al. 
(2018) sampled and sequenced the vaginal microbiome of 
1107 rural Malawi women after pregnancy, and found that 
75.7% (752/994) of the population were dominated by 
Gardnerella vaginalis rather than by Lactobacillus spp., and 
although L. iners increased with time after delivery, G. 
vaginalis still dominated for an extended period. In Doyle’s 
study, both the pregnancy delivery mode and ethnicity also 
appeared to influence the composition of vaginal microbiome, 
though all hosts were healthy. 

Previous research has revealed that the biodiversity of 
vaginal microbial communities varies with health status and 
lifestyle of the host. Nevertheless, existing studies have not 
addressed diversity scaling (changes) across individuals in a 
cohort or population. Theoretically, microbiome diversity 
distribution across individual subjects (i. e., space) is 
traditionally a focus of microbial biogeography. Practically 
speaking, understanding the biogeography of the human 
microbiome can reveal critical information on its 
characteristics in a cohort setting, which can, in turn, 
significantly influence studies on the etiology and 
epidemiology of human microbiome-associated diseases such 
as inflammatory bowel disease, obesity, and BV. To effectively 
assess the spatial scaling of human vaginal microbial 
diversity, we applied the DAR model, which is a recent 
extension of the classic SAR in biogeography and 
conservation biology (Bell et al., 2005; Horner-Devin et al., 
2004; MacArthur & Wilson, 1967; Noguez et al., 2005; Peay et 
al., 2007; Triantis et al., 2012; Varbird et al., 2017; Whittaker 
& Triantis, 2012). SAR is one of the oldest described 
ecological laws or patterns, whereby species richness 
increases with increasing sampling area, and can be traced 
back to the 19th century (Watson, 1835). It is still considered 
one of the most important principles in conservation biology 
and biogeography. The extensions from SAR to DAR 
introduced a _ several important advances including: 
(1) Expanding species richness (number of species) to 
general diversity measures in Hill numbers (Chao et al., 2012, 
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2014; Hill, 1973), thus making it possible to not only assess 
the scaling of species richness (numbers), but also scaling of 
general diversity (e. g., change in community evenness or 
dominance). Therefore, the classic SAR is a special case of 
the more general DAR; (2) The DAR, PDO, MAD, and local 
regional/global diversity (LRD/LGD) profiles are effective tools 
for the biogeographic mapping of biodiversity over space (Ma, 
2018a, 2018c, 2019). 

In this study, we applied DAR modeling and associated 
biogeographic profiles to investigate the spatial diversity 
scaling of postpartum vaginal microbial communities across 
individuals by reanalyzing the large vaginal microbiome 
dataset originally reported by Doyle et al. (2018). The spatial 
diversity scaling of the vaginal microbiome revealed 
heterogeneity among individuals, which could provide an 
ecological basis for personalized and precise diagnosis and 
treatment of microbiome-associated diseases, including BV. 
The biogeographic profiles of the vaginal microbiome also 
provide tools for explaining the DSR hypothesis for BV 
etiology from multiple dimensions (Ma & Ellison, 2018, 2019). 

The vaginal microbial dataset (Doyle et al., 2018) 
reanalyzed in this study consisted of 1 158 vaginal microbiome 
samples collected from 1107 rural Malawi women post- 
delivery. Most samples were collected within the first 20 d of 
delivery, though some were sampled 5-583 d post-delivery. 
The V5-V7 hypervariable regions of the 16S rRNA genes 
were amplified and sequenced under the MiSeq Illumina 
platform. After quality control, the sequences were clustered 
into 14 354 operational taxonomic units (OTUs) using QIIME 
2.8.6. Samples with less than 2 000 reads were removed, as 
were OTUs with less than 1 000 reads. After prescreening, 
1076 samples and 466 OTUs remained for DAR analysis. In 
DAR analysis, the number of each OTU read is equivalent to 
the population abundance of a species in macro-ecology, or 
OTU abundance in diversity analysis. More detailed 
information on the dataset can be found in Doyle et al. (2018). 

The Hill numbers (Hill, 1973) were reintroduced to ecology 
by Jost (2007) and Chao et al. (2012, 2014), and possess 
certain critical advantages over traditional diversity indexes. 
The Hill numbers for measuring alpha diversity are as follows: 


s\MtHa) 
‘D= | S pr] (1) 
i=1 


When q=1, the Hill number is undefined, but its limit exists in 
the following form: 


S 
'D= lim °D = p|- 2 pios(p,)| (2) 


where, D is the diversity in Hill numbers, q (=0, 1, 2,...) is the 
order number of diversity, S is the number of species (or 
OTUs), and p; is the relative abundance of species i. The 
diversity order (q) sets the sensitivity of the Hill numbers to the 
relative frequencies of species abundances. When q=0, °D is 
equal to the number of species or species richness (S). When 
q=1, 1D is the number of typical or common species in the 
community and is equal to the exponential of Shannon 
entropy. When q=2, 2D is more sensitive to species with high 
abundance, and is equal to the inverse of the Simpson index. 


Generally, 1D is the diversity of a community with x="D equally 
abundant species. 

Beta-diversity can be defined with the multiplicative 
partitioning of Hill numbers (Chao et al., 2012, 2014; Ellison, 
2010; Gotelli & Chao, 2013; Jost, 2007), as follows: 

°D; = °D,/ °D, (3) 
where, 1D, and D, are the alpha and gamma diversities in 
terms of Hill numbers, respectively. As 7D, is equivalent to the 
alpha diversity of the meta-community, it has the same 
definition as alpha-diversity (Eqn. (1)). Chao et al. (2012, 
2014) defined a series of Hill numbers corresponding to 
different diversity orders (q) as the diversity profile. In this 
study, the diversity or Hill numbers were computed until the 
third order, g=3. 

According to Ma (2018a), we used the power law (PL) DAR 
model and power-law with exponential cutoff (PLEC) model as 
the DAR models for the human vaginal microbiome. The PL 
model is: 

"D = cA’ (4) 
where, 1D is diversity measured in Hill numbers of the q-th 
order, A is the area (number of individuals), and c & z are the 
PL parameters. 

The PLEC model is: 

1D = cA’exp(dA), (5) 
where, d is a third parameter that is usually less than zero in 
DAR modeling, and exp(dA) is then the exponential decay 
item that eventually overwhelms the power law behavior when 
Ais sufficiently large. 

To simplify parameter estimation, we transformed non-linear 
Equations (4) and (5) into log-linear regression equations: 

In(D ) = In(c) + zin(A) (6) 
In(D) = In(c) + zin(A)+ dA (7) 

In Eqn. (6), z is the slope of the log-linear transformed PL 
model, which is equivalent to its interpretation in the traditional 
SAR—ratio of diversity accrual rate to area increase rate. 
Parameter c of the PL model can be viewed as the number of 
species equivalent to diversity in the first unit of area to 
accrue. Thus, the accrual order of area unit may influence 
parameter c. To deal with this technical issue, the units 
(individuals/samples) to be accumulated were randomly 
permutated each time the DAR model was built. For each 
dataset, we repeatedly applied DAR modeling 100 times by 
randomly re-ordering all samples in the dataset. For the 
detailed computational procedure, please refer to Ma (2018a). 

Similar to the diversity profile concept of Chao et al. (2012, 
2014), which is a series of Hill numbers corresponding to 
different diversity orders (q), Ma (2018a) and Ma & Li (2018) 
proposed four DAR-based profiles, including the DAR, PDO, 
MAD, and LRD/LGD profiles. These four profiles can be 
quantitatively characterized by parameters from the PL/PLEC 
DAR models and can be used to sketch out biogeography 
maps of the human microbiome or other ecological 
communities. 

The DAR profile was defined as a series of z-values 
(scaling parameter) of the PL-DAR model (Eqns. 4 & 6), i.e., a 
series of z-values corresponding to different diversity orders 
(q) or z-q trends. 


The PDO profile was defined as: 
g = 2-2 (8) 
where, z is the scaling parameter of the PL-DAR model, i.e., 
the PDO profile is a series of g (q) values corresponding to 
different diversity orders (q), computed with Eqn. (8). 
The MAD profile was defined as a series of MAD or "Dmax 
values, corresponding to different diversity orders (q): 
"Dmax = e[-ž) exp(—-z ) = cA, exp(-z ) (9) 


max 


where, A,,,, =—z/d is the number of individuals (samples) 
needed to reach the MAD, and c and z are parameters of the 
PLEC-DAR model (Eqns. (5) & (7)). 

The RIP profile was defined as a series of RIP values 
corresponding to different diversity orders (q), as specified by 
the following equation: 

RIP = “cl °D (10) 
where, c is a parameter of the PL-DAR model and D is the 
diversity in Hill numbers estimated with the PLEC-DAR model 
(Eqns. (5) & (7)). Based on the above RIP definition, a RIP 
profile can be defined for a population (cohort) of any size. In 
practice, using Dmax for 7D is more convenient, i.e.: 

TRIP = “cl “Dmax . (11) 

The RIP parameter assesses the average level of an 
individual to represent a population (or cohort) from which the 
individual is a member. The RIP profile is also known as the 
LRD  (local-to-regional diversity) or LGD (local-to-global 
diversity) profile in other ecological systems beyond the 
human microbiome (Ma & Li, 2018; Ma, 2019). 

We built two DAR models for the vaginal microbiome, 
including the PL and PLEC models for alpha-diversity and 
beta-diversity scaling, respectively. The results are listed in 
Tables 1 and 2, including the diversity order (q) of Hill 
numbers, mean model parameters (z, In(c), d, g, Dmax) and 
their standard errors, and measures (correlation coefficient R 
& P-value) for goodness-of-fitting. N represented the number 
of successful fittings out of 100 p re-samplings, as explained 
previously. Re-sampling was performed to deal with the 
possible influence of the order of diversity accrual (i.e., order 
in which the samples were accrued for building the DAR 
model) on model parameter c. Except for two cases of alpha- 
DAR modeling at diversity order g=3, the fittings to the DAR 
models were successful in all 100 re-samplings. Even in the 
two exceptions, the success rates were 97% and 99%, 
respectively. Therefore, the DAR models were considered 
suitable for vaginal microbiome assessment, as also evident 
by the R (linear correlation coefficient) and associated p 
values, which indicated the goodness-of-fit of the DAR models. 

Based on Table 1, we found the following in regard to alpha- 
DAR scaling: 

(1) As one of the most important parameters from the PL- 
DAR model, the scaling parameter (z) at different diversity 
orders (q) was z(q)=(0.807(0), 0.171(1), 0.110(2), 0.095(3)), 
where z(q) represents the DAR profile according to previous 
definition. The DAR profile characterizes the diversity scaling 
across individuals (over space) comprehensively. Results also 
showed that the scaling level differed at different orders. For 
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Table 1 Alpha-DAR models computed with 100 re-samplings for the vaginal microbiome 


Diversity order Power law (PL) 


PL with exponential cutoff (PLEC) 














and statistics z In(c) R P-value g N z d In(c) R P-value N A, Drax 
Mean 0.807 3.902 0.998 0.000 0.250 0.778 0.0001 4.019 0.999 0.000 
SE 0.031 0.200 0.001 0.000 0.038 0.062 0.0001 0.313 0.001 0.000 

q=0 100 100 N/A N/A 
Min 0.738 3.474 0.995 0.000 0.160 0.642 —0.0002 3.201 0.995 0.000 
Max 0.880 4.356 1.000 0.000 0.332 0.948 0.0004 4.683 1.000 0.000 
Mean 0.171 3.345 0.860 0.000 0.874 0.291 -0.0004 2.861 0.930 0.000 
SE 0.039 0.257 0.073 0.000 0.031 0.074 0.0002 0.374 0.044 0.000 

q=1 100 100 667 86.8 
Min 0.085 2.669 0.582 0.000 0.789 0.083 —0.0008 1.820 0.782 0.000 
Max 0.276 3.951 0.963 0.000 0.940 0.488 0.0001 3.831 0.990 0.000 
Mean 0.110 2.808 0.681 0.000 0.920 0.229 —0.0004 2.328 0.829 0.000 
SE 0.047 0.312 0.173 0.000 0.035 0.080 0.0002 0.411 0.107 0.000 

q=2 100 100 530 34.4 
Min 0.013 1.901 0.155 0.000 0.814 0.023 —0.0008 1.162 0.443 0.000 
Max 0.246 3.457 0.933 0.000 0.991 0.448 0.0002 3.225 0.974 0.000 
Mean 0.095 2.565 0.589 0.013 0.931 0.209 —0.0004 2.109 0.770 0.001 
SE 0.054 0.362 0.231 0.084 0.040 0.088 0.0002 0.447 0.157 0.014 

q=3 97 99 507 24.5 
Min —0.002 1.508 0.013 0.000 0.811 —0.004 -0.0009 0.890 0.061 0.000 
Max 0.250 3.269 0.932 0.667 1.001 0.431 0.0003 3.112 0.971 0.138 


z, c & d: The parameters of PL- and/or PLEC-DAR models, in which z of PL-DAR model is the scaling parameter, and c can be viewed as the 
number of species equivalent to diversity in the first unit of area to accrue. g: The pair-wise diversity overlap. R: The correlation coefficient to judge 
the goodness-of-fitting of PL or PLEC models. P-value: The parameter to judge the success or failure of the fitting of PL or PLEC models (Ps0.05 


indicates successful model fitting). N: The number of successful fittings out of 100 re-samplings. A 


: The number of accrued individuals 


max" 


corresponding to the maximal accrual diversity. Dmax: The maximal accrual diversity. SE: Standard error. Min: Minimum. Max: Maximum. N/A: Not 


available. 


example, scaling at diversity order g=0, which is equivalent to 
the classic SAR law, was faster than that at q=1, 2, or 3, as 
C 


0. 
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indicated by the monotonically decreasing z-value (see Figure 
1A for alpha-DAR profile). 
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2 
Diversity order (q) 


Diversity order (q) 


Figure 1 DAR profile (z-q) and PDO profile (g-q) of the vaginal microbiome 


A: Alpha-diversity scaling; B: Beta-diversity scaling. 


(2) The PDO profile was g(q)=(0.250(0), 0.874(1), 0.920(2), 
0.931(3)). The PDO profile, which characterizes the overlap or 
similarity between pair-wise individuals, showed the opposite 
trend as the DAR profile, i.e., a monotonically increasing trend 
(see Figure 1A for alpha-PDO profile). 

(3) The MAD profile characterizes the theoretically maximal 
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accumulation of diversity across individuals. Here, regarding 
the MAD profile, the PLEC model failed to produce D,,,, at 
diversity order q=0 because d>0, for which a maximum does 
not exist. For the diversity orders q=1, 2, 3, the PLEC model 
for alpha-diversity successfully generated D,,,,, i-€., Dma(@) = 
(86.8(1), 34.4(2), 24.5(3)). 


Table 2 Beta-DAR models computed with 100 re-samplings for the vaginal microbiome 


Diversity order Power law (PL) 


PL with exponential cutoff (PLEC) 














and statistics z In(c) R P-value g N z d In(c) R P-value N — Apax Drax 
Mean 0.805 -0.070 1.000 0.000 0.253 0.770 0.00012 0.070 1.000 0.000 
SE 0.010 0.065 0.000 0.000 0.012 0.016 0.00004 0.078 0.000 0.000 

q=0 100 100 N/A N/A 
Min 0.782 -0.217 0.999 0.000 0.228 0.722 -0.00002 -0.131 1.000 0.000 
Max 0.826 0.070 1.000 0.000 0.280 0.812 0.00022 0.313 1.000 0.000 
Mean 0.176 1.591 0.884 0.000 0.870 0.298 -0.00043 1.092 0.952 0.000 
SE 0.031 0.208 0.059 0.000 0.024 0.041 0.00011 0.219 0.034 0.000 

q=1 100 100 687 15.5 
Min 0.104 1.047 0.686 0.000 0.807 0.215 -0.00076 0.449 0.829 0.000 
Max 0.255 2.093 0.974 0.000 0.925 0.435 -0.00022 1.570 0.994 0.000 
Mean 0.146 1.890 0.734 0.000 0.893 0.291 -0.00052 1.295 0.867 0.000 
SE 0.050 0.333 0.178 0.000 0.038 0.076 0.00024 0.368 0.111 0.000 

q=2 100 100 563 17.3 
Min 0.017 1.026 0.106 0.000 0.791 0.102 -0.00124 0.227 0.335 0.000 
Max 0.274 2.730 0.978 0.001 0.988 0.525 -0.00002 2.356 0.987 0.000 
Mean 0.166 1.963 0.711 0.000 0.877 0.330 -0.00058 1.293 0.850 0.000 
SE 0.062 0421 0.200 0.000 0.049 0.103 0.00035 0.468 0.129 0.000 

q=3 100 100 566 21.3 
Min 0.021 0.798 0.108 0.000 0.735 0.069 -0.00149 0.115 0.203 0.000 
Max 0.339 2917 0.972 0.000 0.986 0.591 0.00016 2663 0.980 0.000 


z, c & d: The parameters of PL- and/or PLEC-DAR models, in which z of PL-DAR model is the scaling parameter, and c can be viewed as the 
number of species equivalent to diversity in the first unit of area to accrue. g: The pair-wise diversity overlap. R: The correlation coefficient to judge 
the goodness-of-fitting of PL or PLEC models. P-value: The parameter to judge the success or failure of the fitting of PL or PLEC models (Ps0.05 


indicates successful model fitting). N: The number of successful fittings out of 100 re-samplings. A 


The number of accrued individuals 


max" 


corresponding to the maximal accrual diversity. Dmax: The maximal accrual diversity. SE: Standard error. Min: Minimum. Max: Maximum. N/A: Not 


available. 


From Table 2, we found the following in regard to beta-DAR 
scaling: 

(1) As one of the most important parameters of the PL-DAR 
model, the beta-diversity scaling parameter (z) at different 
diversity orders (q) was z(q) = (0.805(0), 0.176(1), 0.146(2), 
0.166(3)), where z(q) represents the DAR profile according to 
previous definition and characterizes diversity scaling across 
individuals (over space) comprehensively (see Figure 1B for 
beta-DAR profile). Comparison between the beta-DAR and 
alpha-DAR profiles revealed an interesting phenomenon: i.e., 
the alpha-DAR profile monotonically decreased with q, 
whereas the beta-DAR profile was valley-shaped. This 
suggests that, at a lower diversity order (q), the alpha-DAR 
and beta-DAR scaling parameters (z) were rather close to 
each other, but the difference was enlarged at higher diversity 
orders (q). 

(2) The beta-PDO profile was g(q) = (0.253(0), 0.870(1), 
0.893(2), 0.877(3)) for beta-diversity scaling. Here, the PDO 
profile, which characterizes the overlap or similarity between 
pair-wise individuals, showed the opposite trend to the DAR 
profile, i.e., a bell-shaped trend (see Figure 1B for beta-PDO 
profile). 

(3) The MAD profile characterizes the theoretical maximal 
accumulation of diversity across individuals. Here, regarding 
the beta-MAD profile, the PLEC model failed to produce Dmax 
at diversity order g=0 because d>0, for which a maximum 


does not exist. For diversity orders q=1, 2, 3, the PLEC model 
for beta-diversity successfully generated D,,,,, i. €., Dmax(Q) = 
(15.5(1), 17.3(2), 21.3(3)). 

Table 3 shows the RIP values for both alpha-DAR and beta- 
DAR of the vaginal microbiome. At diversity order g=0, the 
estimation of "Das failed, and RIP for q=0 could not be 
estimated. For q=1, 2, 3, RIP was successfully estimated for 
alpha- and beta-diversity, respectively. Here, RIP 
characterized the relationship between individual- and 
population-level diversity. For example, at diversity order q=1, 
alpha-RIP=0.327 and beta-RIP=0.316, indicating that an 
average individual represented approximately 33% and 32% 
of population alpha- and beta-diversity, respectively. 

In the current study, we investigated the diversity (including 
alpha- and beta-diversity) scaling of the human vaginal 
microbiome across individuals by re-analyzing a big dataset 


Table 3 Ratio of individual-level to population-level diversity 
(RIP) of the vaginal microbiome 


Diversity order RIP for alpha-DAR RIP for beta-DAR 





q=0 N/A N/A 

q=1 0.327 0.316 
q=2 0.482 0.383 
q=3 0.530 0.335 


N/A: Not available. 
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originally published by Doyle et al. (2018). Compared with the 
microbial SAR range reported in existing literature for other 
microbes, such as Green & Bohannan’s (2006) range 
between 0.019-—0.470, the scaling parameter (z) estimated in 
our study, i.e., a/pha-z=0.807, beta-z=0.805, appears to be out 
of the known range, at nearly twice that reported for SAR 
values for other microbes. Three possibilities exist for the 
significant difference: (1) The use of revolutionary 
metagenomic sequencing technology, which allows for 
detection of more microbial species and consequently large 
scaling parameter; (2) The human vaginal microbiome has 
higher heterogeneity across individuals, which could be 
validated by future biomedical studies; and (3) The 
postpartum nature of the vaginal microbiome samples 
analyzed in this study. We could not exclude these 
possibilities at present due to insufficient available data for 


comparative research. Indeed, previous studies have 
classified human vaginal microbiomes into five main 
community-state types (CSTs), in which CST |, Il, Ill, and V 


are dominated by Lactobacillus spp., and CST IV is composed 
of facultative or strictly anaerobic bacteria (Gajer et al., 2012; 
Ravel et al., 2011), many of which are BV-related. The vaginal 
microbial communities of postpartum women in rural Malawi 
studied by Doyle et al. (2018) and reanalyzed here were 
mostly Lactobacillus-deficient microbiomes, which could be 
grouped as CST IV, although all these women were healthy. 
Therefore, the classification of CSTs may be more complex 
than initially conceived. Consequently, our DAR analysis 
based on Doyle et al. (2018) may be limited by the datasets of 
postpartum women, and the DAR parameters of the vaginal 
microbiomes of other CST women are likely different from the 
results reported here. Further studies should be performed to 
clarify this important issue. 

The major findings in this study can be summarized using 
four profiles: i. e., DAR profile, characterizing the change 
(scaling) in diversity heterogeneity across individuals; PDO 
profile, characterizing the pair-wise similarity (overlap) 
between individuals; MAD profile, characterizing the maximal 
accrual diversity in a population; and RIP profile, 
characterizing the ratio of individual-level diversity to 
population-level diversity. Theoretically, the four profiles can 
together summarize the essential characteristics of the spatial 
distribution of vaginal microbial diversity and offer effective 
tools to sketch out the biogeographic maps of the human 
vaginal microbiome. Practically, they are essentially 
quantitative metrics of diversity heterogeneity across 
individuals from different dimensions (diversity scaling, pair- 
wise similarity in diversity, maximal accrual diversity, ratio of 
individual to population diversity). These multidimensional 
metrics could provide more comprehensive tools for 
understanding the implications of vaginal microbial diversity to 
women’s health, including the DSR hypothesis for BV etiology 
(Ma et al., 2012, Ma & Ellison 2018, 2019 Sobel, 1999). In 
addition, the quantitative models of the four profiles obtained 
here could be harnessed to assess and predict microbiome 
diversity changes at the population scale and are of potential 
significance for evaluating women’s health associated with 
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vaginal microbiomes. 
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